function [xplates,gz]=calcgrav(filenumber,Nplates,bouguer)
%bouguer ==1 : calculates bouguer anomaly
%bouguer ==0 : calculates free air anomaly

count=filenumber;

N=readdata(count,1);
S=readdata(count,2);
T=readdata(count,3);
%     multigridlevels=readdata(count,4);
%     bound=readdata(count,5);
time=readdata(count,6);
%     dtnom=readdata(count,7);
dx=readdata(count,8);
dz=readdata(count,9);
mx=readdata(count,10);
mz=readdata(count,11);
%     mcs=readdata(count,12);
%     mct=readdata(count,13);
mtype=readdata(count,14);
mTemp=readdata(count,15);
%mstressxx=readdata(count,16);
%mstressxz=readdata(count,17);
%mstrainxx=readdata(count,18);
%mstrainxz=readdata(count,19);
%mfinitestrain=readdata(count,20);
%     vx=readdata(count,21);
%     vz=readdata(count,22);
%P=readdata(count,23);
Nsurfaces=readdata(count,24);
surfaces=readdata(count,25);

if bouguer ==1 sealevel=max(surfaces(:,Nsurfaces)); end
if bouguer ==0 sealevel=20e3; end

%sealevel=surfaces(1,3); %!!!!!!!!!!!!!!!!!!!!!!
%sealevel=max(surfaces(:,Nsurfaces));
%sealevel=min(surfaces(:,Nsurfaces));
%sealevel=20e3;

x=zeros(1,S); z=zeros(1,T);
for s=1:S-1 x(s+1)=x(s)+dx(s); end
for t=1:T-1 z(t+1)=z(t)+dz(t); end


mrho=0*mtype;
ind=find(mtype==1); mrho(ind)=2800-3.28e-5*mTemp(ind)*2800;
ind=find(mtype==2); mrho(ind)=2800-3.28e-5*mTemp(ind)*2800;
ind=find(mtype==3 | mtype==4); mrho(ind)=3300*exp(-2.77e-5*(mTemp(ind)-273)-0.5*0.97e-8*(mTemp(ind).*mTemp(ind)-273*273)-0.32*(1./mTemp(ind)-1/273));
ind=find(mtype==5); mrho(ind)=2800-3.28e-5*mTemp(ind)*2800; mrho(ind)=2300-3.28e-5*mTemp(ind)*2300;
ind=find(mtype==6); mrho(ind)=2800-3.28e-5*mTemp(ind)*2800;
%clear('mtemp','mtype');

[grid,xgrid,zgrid]=bilingrid(mrho,mx,mz,dx,dz,1,1); zgrid=zgrid-sealevel; [Zgrid,Xgrid] = meshgrid(zgrid,xgrid);

Splates=1024; zmaxplates=500e3; dxplates=max(x)/Splates; xplates=linspace(dxplates/2,max(x)-dxplates/2,Splates);
H=ones(1,Nplates); H=linspace(5,50,Nplates); H=H/sum(H)*zmaxplates; rho=zeros(length(H),Splates);
zpos=H(1)/2; Xplates=0*rho; Zplates=0*rho;
for h=1:length(H)
    Xplates(h,:)=xplates;
    Zplates(h,:)=zpos*ones(1,Splates);
    if h<length(H)
        zpos=zpos+(H(h)+H(h+1))*0.5;
    end
end
rho=interp2(Xgrid',Zgrid',grid',Xplates,Zplates);
rhomean=mean(rho'); deltarho=0*rho;
for h=1:length(H)
    deltarho(h,:)=rho(h,:)-rhomean(h);
end

gz=grav1D(H,deltarho,sum(dx),zmaxplates);



end